subroutine SetInitConditions()
use Data_cfoam
use Data_Mesh
use Acoustics, only: AcousticSurf, nNoiseSignals, NoiseSignal
use DataTime, only: nStepInitial, ttime, dt,dt2, AsynchronousTime, TimeCell,dtCell,dt2Cell,dtFace, &
                    TimeStepMin,TimeStepMax, TimeLocal, CurrentTimeCell, dtFaceCorr, NextTimeCell, &
                    NextTimeLocal, nNextTimeCell,nNextTimeLocal, nNextTimeCell1, &
                    nTimeStepDiff, TimeStepDiff,TimeStepMin, dnTimeCell
use DataCells
use DataFaces
use DataPerfectGas
use Data_Transducer, only: nTransducers, Transducer
use DataViscosity, only: CalcViscosityOn
use InletBoundaryConditions, only: InitCondInletBCon
implicit none
    include 'mpif.h'

    integer :: i,j,k, jj, iCell,iFace
    real(8) :: xg(3)
    real(8) :: xg1  ! node's coord for postprocessing
    real(8) :: xg5  ! node's coord for postprocessing
    real(8) :: rr1, xt(3,8)
    integer :: nCellsOutSum !, nUnusedSectionsSum ! for parallel
    character(3) :: DirString
    real(8) :: cZ(3) ! coord of cell centers

    integer :: iNS, iS,iL, nFL, iTrans
    character(280) :: nameFile
    
    real(8) :: auxReadReal

    real(8) :: c,r(3),Q,dQdT,P0,U0(3),Rho0, TimeX, rMod
    real(8) :: x,y,z, a,b

    real(8) :: L,eps,alfa, x0,y0,tau,Tetta, coef1,coef2, rr, Umean,Px, sigma,Pi,U00
    
    if(rank.eq.0)then
        write(*,*) 'Initialization ...'
    end if
    
    nStepInitial = 1

	allocate (p(nCells+nFaces-nInternalFaces))
	allocate (T(nCells+nFaces-nInternalFaces))
	allocate (U(3,nCells+nFaces-nInternalFaces))
	allocate (Ps(nFaces))
	allocate (Ts(nFaces))
	allocate (Us(3,nFaces))
	allocate (pnew(nCells+nFaces-nInternalFaces))
	allocate (Tnew(nCells+nFaces-nInternalFaces))
	allocate (Unew(3,nCells+nFaces-nInternalFaces))
	allocate (Psnew(nFaces))
	allocate (Tsnew(nFaces))
	allocate (Usnew(3,nFaces))
	allocate (gradU(3,3,nPoints)) !nFaces))
!	allocate  (meshCellOut(nCells))
	allocate (meshCellOut(nCells),areaOut(nCells),xOut(3,nCells),normalOut(3,nCells)) ! Grigoriy's new acoustics postprocessor
	allocate (rho(nCells+nFaces-nInternalFaces))           !
	allocate (E(nCells+nFaces-nInternalFaces))             !
	allocate (Rhos(nFaces))                                !
	allocate (Es(nFaces))                                  !
	allocate (rhonew(nCells+nFaces-nInternalFaces))        !
	allocate (Enew(nCells+nFaces-nInternalFaces))          !
	allocate (Rhosnew(nFaces))                             !
	allocate (Esnew(nFaces))                               !
	allocate(VarMem(nFaces))


    ! initialization
!	nu=0.0d0
!	gam=1.4d0
	cv=717.65d0 ! [0 2 -2 -1 0 0 0] 
	U=0.0d0
	nCellsOut=0
	nUnusedSections=0
!        ttime=0.0d0
!	call interp(p1,np1,ttime,uin)

    call CalcGammas()

    ! for postprocessing
    if(AcousticsOnOff) then
    
        if(AcousticSurf)then
            call AcousticSurfInit()
            call NoiseSugnalInit()
!            call ReadContourData('ContourData.dat')
!        end if
        else
	        xg1=1e10
	        xg5=1e10
	        do i=1,nCells
	            xg(1:3)=0.0d0
	            rr1=-1d0
	            do j=1,8
	                xg(1:3)=xg(1:3)+meshPoints(1:3,meshCellPoints(j,i))/8.0d0
	                if (dabs(meshPoints(2,meshCellPoints(j,i))).lt.1d-5) rr1=1d0
	            end do
        	    
        !        include 'meshCellOut.H' ! Grigoriy's acoustics postprocessor
                do jj=1,8
	                xt(1:3,jj)=meshPoints(1:3,meshCellPoints(jj,i))
	            end do			
	            call meshCellOut_defn(xt,meshCellOut,nCellsOut,i,&
		                areaOut(nCellsOut+1),xOut(1:3,nCellsOut+1),normalOut(1:3,nCellsOut+1), &
		                nUnusedSections)

        ! instead of this:
        !        rr=dsqrt(xg(2)**2+xg(3)**2)
        !	    
        !	    rr1=0.158+0.158*14.0/180.0*3.1415927*xg(1)
        !        rr2=0.168+0.168*14.0/180.0*3.1415927*xg(1)
        !        if (xg(1).ge.0.0.and.rr.ge.rr1.and.rr.le.rr2) then
        !            nCellsOut=nCellsOut+1
        !            meshCellOut(nCellsOut)=i
        !        end if
        ! ------------------------------------------------------------------------

                !!!        if (xg(2).ge.0.0.and.rr1.gt.0.0d0) then
                !!!        nCellsOut=nCellsOut+1
                !!!        meshCellOut(nCellsOut)=i
                !!!        end if
                                                                
                !	if (dabs(xg(1)-0.05d0).lt.xg1 &
                !.and.dabs(xg(2)).lt.5e-4.and.xg(3).ge.1.0d-5.and.xg(3).lt.0.1d0) &
                !	xg1=dabs(xg(1)-0.05d0)

                !	if (dabs(xg(1)-0.25d0).lt.xg5 &
                !.and.dabs(xg(2)).lt.5e-4.and.xg(3).ge.1.0d-5.and.xg(3).lt.0.1d0) &
                !	xg5=dabs(xg(1)-0.25d0)
	        end do

            if (nump.gt.1) then
	            nCellsOutSum = 0
	            nUnusedSectionsSum = 0
                call MPI_ALLREDUCE(nCellsOut,nCellsOutSum,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
	            call MPI_Barrier(MPI_COMM_WORLD,ierr)

                call MPI_ALLREDUCE(nUnusedSections,nUnusedSectionsSum,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
	            call MPI_Barrier(MPI_COMM_WORLD,ierr)
            else
   	            nCellsOutSum = nCellsOut
	            nUnusedSectionsSum = nUnusedSections
            end if
            
            if(rank.eq.0)then
                if(nCellsOutSum.le.0)then
                    stop '@@Fatal Error: wrong acoustic surface: nCellsOut==0!'
                else
                    if(rank.eq.0)then
                        write(*,*) 'Surface: nCellsOut = ',nCellsOutSum, ' nUnusedSections = ',nUnusedSectionsSum
                        nUnusedSectionsSum = nUnusedSectionsSum * 100 / nCellsOutSum
                        if(nUnusedSectionsSum.ne.0)then
                            write(*,*) '@Warning: nUnusedSections = ',nUnusedSectionsSum, '%'
                        else
                            write(*,*) 'nUnusedSections = 0%'
                        end if
                    end if
                end if
            end if

	    end if
    end if
    
!	write(*,*) 'rank,xg1,xg5=',rank,xg1,xg5
    ! end for postprocessing

!	call inputScalar(processor(1:ie)//'0/T',nCells,T)
!	call inputScalar(processor(1:ie)//'0/p',nCells,p)
!	call inputVector(processor(1:ie)//'0/U',nCells,U)

!***************************************************************************
! set initial conditions
!***************************************************************************
    if(irestart.le.0) then
        RoFon = PFON/(TFON*cv*(Gam-1.0d0)) 
	    do i=1,nCells
	        xg(1:3)=0.0d0
	        do j=1,8
	            xg(1:3)=xg(1:3)+meshPoints(1:3,meshCellPoints(j,i))/8.0d0
	        end do
    	
    !	if (dabs(xg(2)).lt.5e-4 ) then
    !	nCellsOut=nCellsOut+1
    !	meshCellOut(nCellsOut)=i
    !	end if
    !	call vikhr(xg,U(1:3,i),p(i),rho(i),E(i))
    !    T(i)=p(i)/(rho(i)*cv*(Gam-1.0d0))
	        cZ(1:3) = meshCellCenters(1:3,i) 
	        U(1,i)= UFON ! 0d0 ! cZ(1) !0d0 !
	        U(2,i)= 0d0 ! cZ(2) !0d0 !VFON+omega*cZ(3)
	        U(3,i)= 0d0 ! cZ(3) !0d0 !WFON-omega*cZ(2)
	        p(i)=PFON
	        t(i)=TFON
	        
!            U(1:3,i) = 0.0
!            do k=1,6
!                if(meshFaceOwner(meshCells(k,i)).eq.i)then
!                    U(1,i) = U(1,i) + meshSf(1,meshCells(k,i))
!                    U(2,i) = U(2,i) + meshSf(2,meshCells(k,i))
!                    U(3,i) = U(3,i) + meshSf(3,meshCells(k,i))
!                else
!                    U(1,i) = U(1,i) - meshSf(1,meshCells(k,i))
!                    U(2,i) = U(2,i) - meshSf(2,meshCells(k,i))
!                    U(3,i) = U(3,i) - meshSf(3,meshCells(k,i))
!                end if
!            end do


        !	rr=dsqrt(xg(2)**2+xg(3)**2)
        !	if (rr.le.0.025d0) u(1,i)=UFON
        	
        !!    p(i)=1.0d5
        !    T(i)=287.0d0
        !!	rho(i)=1d0                                                            
    !	    rho(i)=p(i)/(T(i)*cv*(Gam-1.0d0))                                     !!
        !!	T(i)=p(i)/(rho(i)*cv*(Gam-1.0d0))                                     !!

    	    rho(i)=p(i)/(T(i)*cv*(Gam-1.0d0)) 
	        E(i)=p(i)/((gam-1.0d0)*rho(i))+(U(1,i)**2+U(2,i)**2+U(3,i)**2)/2.0d0  !! velocities is zero before it
	        E(i) = T(i)*cv + (U(1,i)**2+U(2,i)**2+U(3,i)**2)/2.0d0
	    end do
    	
	    ! set primary variables on surfaces
        do k=1,nInternalFaces
	        Psnew(k)=p(meshFaceOwner(k))
	        Tsnew(k)=T(meshFaceOwner(k))
	        Esnew(k)=E(meshFaceOwner(k))        !!
	        Rhosnew(k)=rho(meshFaceOwner(k))    !!
	        Usnew(1:3,k)=U(1:3,meshFaceOwner(k))
        end do
        
	    do j=1,nBoundaryFields
	        do i=1,boundaryMeshSize(j)
	            Psnew(i+boundaryMeshStart(j))=p(meshFaceOwner(i+boundaryMeshStart(j)))
	            Tsnew(i+boundaryMeshStart(j))=T(meshFaceOwner(i+boundaryMeshStart(j)))
	            Esnew(i+boundaryMeshStart(j))=E(meshFaceOwner(i+boundaryMeshStart(j)))     !!
	            Rhosnew(i+boundaryMeshStart(j))=rho(meshFaceOwner(i+boundaryMeshStart(j))) !!
	            Usnew(1:3,i+boundaryMeshStart(j))=U(1:3,meshFaceOwner(i+boundaryMeshStart(j)))
	        end do
	    end do

        c = sqrt(1.4*PFon/RoFon)
! Initial conditions from Calc3D
        if(InitCondCalc3D) then
            do iCell = 1,nCells
                x = meshCellCenters(1,iCell)
                c = 340d0
                TimeX = x / c
                r = 0d0
                call CalcQ(r,TimeX,Q,dQdt,p0,U0)

                r(1:3) = meshCellCenters(1:3,iCell)
                rMod = sqrt(sum(r(1:3)*r(1:3)))
                c = sqrt(1.4*PFon/RoFon)
                TimeX = - rMod / c
                call CalcQ3D(r,0d0,Q,dQdt,p0,U0)
                r(1:3) = -meshCellCenters(1:3,iCell)
                call CalcQ3D(r,0d0,Q,dQdt,p0,U0)

                Pnew(iCell) = p0
                Unew(1:3,iCell) = U0(1:3)
                Rhonew(iCell) = RoFon + (P0-PFon)/c**2
        	    Enew(iCell) = P0 / (0.4 * (RoFon + (P0-PFon)/c**2)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
    !            Pnew(iCell) = 0.4d0 * (Enew(iCell) - 0.5d0 * sum(U0(1:3)*U0(1:3))) * Rhonew(iCell)! * 0.4d0

                p(iCell) = Pnew(iCell)
                U(1:3,iCell) = Unew(1:3,iCell)
                E(iCell) = Enew(iCell)
                Rho(iCell) = Rhonew(iCell)
            end do

            do iFace = 1,nFaces
                x = -meshFaceCenters(1,iFace)
                c = 340d0
                TimeX = x / c
                call CalcQ(r,TimeX,Q,dQdt,p0,U0)

                r(1:3) = -meshFaceCenters(1:3,iFace)
                rMod = sqrt(sum(r(1:3)*r(1:3)))
                c = sqrt(1.4*PFon/RoFon)
                TimeX = - rMod / c
                call CalcQ3D(r,0d0,Q,dQdt,p0,U0)
                r(1:3) = meshFaceCenters(1:3,iFace)
                call CalcQ3D(r,0d0,Q,dQdt,p0,U0)

                Psnew(iFace) = p0
                Usnew(1:3,iFace) = U0(1:3)
        	    Esnew(iFace) = P0 / (0.4 * (RoFon + (P0-PFon)/c**2)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
                Rhosnew(iFace) = RoFon + (P0-PFon)/c**2

                ps(iFace) = Psnew(iFace)
                Us(1:3,iFace) = Usnew(1:3,iFace)
                Es(iFace) = Esnew(iFace)
                Rhos(iFace) = Rhosnew(iFace)
            end do
        end if

! Cylindrical waves
        if(InitCondHankelFunc2D) then
            do iCell = 1,nCells
                r(1:3) = -meshCellCenters(1:3,iCell)
!                call CalcQ2D_Hankel(r,0d0,p0,U0,rho0)
                call CalcQ3D_ConvectedMonopole(r,0d0,p0,U0,rho0)

                Pnew(iCell) = p0
                Unew(1:3,iCell) = U0(1:3)
                Rhonew(iCell) = Rho0
        	    Enew(iCell) = P0 / (0.4 * Rho0) + 0.5d0 * sum(U0(1:3)*U0(1:3))
    !            Pnew(iCell) = 0.4d0 * (Enew(iCell) - 0.5d0 * sum(U0(1:3)*U0(1:3))) * Rhonew(iCell)! * 0.4d0

                p(iCell) = Pnew(iCell)
                U(1:3,iCell) = Unew(1:3,iCell)
                E(iCell) = Enew(iCell)
                Rho(iCell) = Rhonew(iCell)
            end do

            do iFace = 1,nFaces
                r(1:3) = meshFaceCenters(1:3,iFace)
!                call CalcQ2D_Hankel(r,0d0,p0,U0,rho0)
                call CalcQ3D_ConvectedMonopole(r,0d0,p0,U0,rho0)

                Psnew(iFace) = p0
                Usnew(1:3,iFace) = U0(1:3)
                Rhosnew(iFace) = Rho0
        	    Esnew(iFace) = P0 / (0.4 * Rho0) + 0.5d0 * sum(U0(1:3)*U0(1:3))

                ps(iFace) = Psnew(iFace)
                Us(1:3,iFace) = Usnew(1:3,iFace)
                Es(iFace) = Esnew(iFace)
                Rhos(iFace) = Rhosnew(iFace)
            end do
        end if
        
! 1D waves
        if(InitCond1Dwaves) then
            do iCell = 1,nCells
                r(1:3) = -meshCellCenters(1:3,iCell)
                call CalcQ(r,0d0,Q,dQdt,p0,U0)
                

                Pnew(iCell) = p0
                Unew(1:3,iCell) = U0(1:3)
                Rhonew(iCell) = RoFon + (P0-PFon)/c**2
        	    Enew(iCell) = P0 / (0.4 * (RoFon + (P0-PFon)/c**2)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
    !            Pnew(iCell) = 0.4d0 * (Enew(iCell) - 0.5d0 * sum(U0(1:3)*U0(1:3))) * Rhonew(iCell)! * 0.4d0

                p(iCell) = Pnew(iCell)
                U(1:3,iCell) = Unew(1:3,iCell)
                E(iCell) = Enew(iCell)
                Rho(iCell) = Rhonew(iCell)
            end do

            do iFace = 1,nFaces
                r(1:3) = meshFaceCenters(1:3,iFace)
                call CalcQ(r,0d0,Q,dQdt,p0,U0)

                Psnew(iFace) = p0
                Usnew(1:3,iFace) = U0(1:3)
                Rhosnew(iFace) = RoFon + (P0-PFon)/c**2
        	    Esnew(iFace) = P0 / (0.4 * (RoFon + (P0-PFon)/c**2)) + 0.5d0 * sum(U0(1:3)*U0(1:3))

                ps(iFace) = Psnew(iFace)
                Us(1:3,iFace) = Usnew(1:3,iFace)
                Es(iFace) = Esnew(iFace)
                Rhos(iFace) = Rhosnew(iFace)
            end do
        end if
        
        if(Benchmark_CAA2)then
            do iCell = 1,nCells
                r(1:3) = -meshCellCenters(1:3,iCell)

                Pnew(iCell) = Pfon + 2d0**( -( (r(1)-4d0)**2+r(2)**2 ) / 0.04d0 )
                Unew(1:3,iCell) = 0d0
                Rhonew(iCell) = RoFon
        	    Enew(iCell) = Pnew(iCell) / (0.4 * Rhonew(iCell)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
        	    Tnew(iCell) = Tfon

                p(iCell) = Pnew(iCell)
                U(1:3,iCell) = Unew(1:3,iCell)
                E(iCell) = Enew(iCell)
                Rho(iCell) = Rhonew(iCell)
                T(iCell) = Tnew(iCell)
            end do

            do iFace = 1,nFaces
                r(1:3) = meshFaceCenters(1:3,iFace)
                call CalcQ(r,0d0,Q,dQdt,p0,U0)

                Psnew(iFace) = Pfon + 2d0**( -( (r(1)-4d0)**2+r(2)**2 ) / 0.04d0 )
                Usnew(1:3,iFace) = 0d0
                Rhosnew(iFace) = RoFon
        	    Esnew(iFace) = Psnew(iFace) / (0.4 * Rhosnew(iFace)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
        	    Tsnew(iFace) = Tfon

                ps(iFace) = Psnew(iFace)
                Us(1:3,iFace) = Usnew(1:3,iFace)
                Es(iFace) = Esnew(iFace)
                Rhos(iFace) = Rhosnew(iFace)
                Ts(iFace) = Tsnew(iFace)
            end do
        end if
        
!         ! Viscosity Checking
!        if(ViscPartOnly)then            
!            L = 0.005
!            Umean = 0.01
!            Px = - 12d0 * nu / L**2 * Umean
!
!            do i = 1,nCells
!                x = -meshCellCenters(1,i)
!                y = -meshCellCenters(2,i) !- 0.5d0
!                z = -meshCellCenters(3,i)
!                P(i) = Pfon + Px * x
!!                U(1,i) = a*(y**2 - 0.25d0)
!                
!                U(1,i) = -0.5d0 / nu * Px * z*(L-z)
!                
!                !U(1,i) = 0.5 * (Pout-Pin) * (y-1d0) * y / nu
!                U(2,i) = 0d0
!                U(3,i) = 0d0
!                T(i) = Tfon
!                Rho(i) = DensityPT(P(i),T(i))
!                E(i) = P(i) / (0.4d0*Rho(i)) + 0.5d0 * U(1,i)**2
!
!                Pnew(i) = P(i)
!                Unew(1:3,i) = U(1:3,i)
!                Enew(i) = E(i)
!                Rhonew(i) = Rho(i)
!                Tnew(i) = T(i)
!            end do
!
!            do i = 1,nFaces
!                x = meshFaceCenters(1,i)
!                y = meshFaceCenters(2,i) !- 0.5d0
!                z = meshFaceCenters(3,i)
!                Ps(i) = Pfon + Px * x
!                Rhos(i) = RoFon
!                Us(1,i) = -0.5d0 / nu * Px * z*(L-z)
!
!                !Us(1,i) = 0.5 * (Pout-Pin) * (y-1d0) * y / nu
!                Us(2,i) = 0d0
!                Us(3,i) = 0d0
!                Ts(i) = Tfon
!                Rhos(i) = DensityPT(Ps(i),Ts(i))
!                Es(i) = Ps(i) / (0.4d0*Rhos(i)) + 0.5d0 * Us(1,i)**2
!
!!!                Ps(i) = 0d0
!!                Rhos(i) = 0d0
!!!                Us(1,i) = 0d0
!!                Us(2,i) = 0d0
!!                Us(3,i) = 0d0
!!                Es(i) = 0d0
!!                Ts(i) = 0d0
!
!                Psnew(i) = Ps(i)
!                Usnew(1:3,i) = Us(1:3,i)
!                Esnew(i) = Es(i)
!                Rhosnew(i) = Rhos(i)
!                Tsnew(i) = Ts(i)
!            end do
!            
!        end if

        ! Viscosity Checking
        if(ViscPartOnlyX)then
            Pin = 1.0001d5
            Pout = 1d5
            RoFon = 1d0

            do i = 1,nCells
                x = -meshCellCenters(1,i) !- 5d0
                y = -meshCellCenters(2,i) !- 0.5d0
                P(i) = Pin + (Pout-Pin) * y * 10d0
                Rho(i) = RoFon
                a = 0.5 * (Pout-Pin) / nu * 10d0
                b = 0d0
                c = -a * 0.25
                U(2,i) = a*x*(x-10d0)
                !U(1,i) = 0.5 * (Pout-Pin) * (y-1d0) * y / nu
                U(1,i) = 0d0
                U(3,i) = 0d0
                E(i) = P(i) / (1.4d0*RoFon) + 0.5d0 * U(1,i)**2
                T(i) = 0d0

                Pnew(i) = P(i)
                Unew(1:3,i) = U(1:3,i)
                Enew(i) = E(i)
                Rhonew(i) = Rho(i)
                Tnew(i) = T(i)
            end do

            do i = 1,nFaces
                x = meshFaceCenters(1,i) !- 5d0
                y = meshFaceCenters(2,i) !- 0.5d0
                Ps(i) = Pin + (Pout-Pin) * y * 10d0
                Rhos(i) = RoFon
                a = 0.5 * (Pout-Pin) / nu * 10d0
                b = 0d0
                c = -a * 0.25
                Us(2,i) = a*x*(x-10d0)
                !Us(1,i) = 0.5 * (Pout-Pin) * (y-1d0) * y / nu
                Us(1,i) = 0d0
                Us(3,i) = 0d0
                Es(i) = Ps(i) / (1.4d0*RoFon) + 0.5d0 * Us(1,i)**2

!                Ps(i) = 0d0
                Rhos(i) = 1d0
                Es(i) = 0d0
                Ts(i) = 0d0

                Psnew(i) = Ps(i)
                Usnew(1:3,i) = Us(1:3,i)
                Esnew(i) = Es(i)
                Rhosnew(i) = Rhos(i)
                Tsnew(i) = Ts(i)
            end do
        end if

        ! Isentropic Vortex
        if(IsentropicVortex)then
            L = 5d0
            eps = 0.3d0
            alfa = 0.204d0
            x0 = 50d0
            y0 = 50d0
            Pfon = 1d0
            RoFon = 1d0

            do i = 1,nCells
                x = -meshCellCenters(1,i)
                y = -meshCellCenters(2,i) !- 0.5d0

                rr = sqrt((x-x0)**2 + (y-y0)**2)
                tau = rr/L
                Tetta = atan2(y-y0,x-x0)
                
                coef1 = (gam-1d0)*0.25d0/(alfa*gam) * eps**2 * exp(2d0*alfa*(1d0-tau**2))
                coef2 = eps*tau * exp(alfa*(1d0-tau**2))
                
                Rho(i) = RoFon * (1d0 - coef1)**(1d0/(gam-1d0))
                U(1,i) = coef2 * sin(Tetta)
                U(2,i) = -coef2 * cos(Tetta)
                U(3,i) = 0d0
                P(i) = Pfon * (1d0 - coef1)**(gam/(gam-1d0))
                E(i) = P(i) / (0.4d0*Rho(i)) + 0.5d0 * (U(1,i)**2 + U(2,i)**2)
                T(i) = TemperaturePRho(P(i),Rho(i))

                Pnew(i) = P(i)
                Unew(1:3,i) = U(1:3,i)
                Enew(i) = E(i)
                Rhonew(i) = Rho(i)
                Tnew(i) = T(i)
            end do

            do i = 1,nFaces
                x = meshFaceCenters(1,i)
                y = meshFaceCenters(2,i) !- 0.5d0

                rr = sqrt((x-x0)**2 + (y-y0)**2)
                tau = rr/L
                Tetta = atan2(y-y0,x-x0)
                
                coef1 = (gam-1d0)*0.25d0/(alfa*gam) * eps**2 * exp(2d0*alfa*(1d0-tau**2))
                coef2 = eps*tau * exp(alfa*(1d0-tau**2))
                
                Rhos(i) = RoFon * (1d0 - coef1)**(1d0/(gam-1d0))
                Us(1,i) = coef2 * sin(Tetta)
                Us(2,i) = -coef2 * cos(Tetta)
                Us(3,i) = 0d0
                Ps(i) = Pfon * (1d0 - coef1)**(gam/(gam-1d0))
                Es(i) = Ps(i) / (0.4d0*Rhos(i)) + 0.5d0 * (Us(1,i)**2 + Us(2,i)**2)
                Ts(i) = TemperaturePRho(Ps(i),Rhos(i))

                Psnew(i) = Ps(i)
                Usnew(1:3,i) = Us(1:3,i)
                Esnew(i) = Es(i)
                Rhosnew(i) = Rhos(i)
                Tsnew(i) = Ts(i)
            end do

        end if
        
        if(GaussVelocityProfile)then
            sigma = 0.025d0
            Pi = 3.1415926535897932384626433832795
            
            do i = 1,nCells
                x = -meshCellCenters(3,i)
                y = -meshCellCenters(2,i) !- 0.5d0
                U00 = 1.0d2

                rr = sqrt((x-x0)**2 + (y-y0)**2)
                if(rr.lt.0.03)then
                    U(1,i) = U00
                else
                    U(1,i) = exp(-0.5d0*((rr-0.03)/sigma)**2) * U00
                end if
                
                Rho(i) = RoFon
                U(2,i) = 0d0
                U(3,i) = 0d0
                P(i) = Pfon
                E(i) = P(i) / (0.4d0*Rho(i)) + 0.5d0 * (U(1,i)**2 + U(2,i)**2)
                T(i) = TemperaturePRho(P(i),Rho(i))

                Pnew(i) = P(i)
                Unew(1:3,i) = U(1:3,i)
                Enew(i) = E(i)
                Rhonew(i) = Rho(i)
                Tnew(i) = T(i)
            end do
            
            do i = 1,nFaces
                x = meshFaceCenters(3,i)
                y = meshFaceCenters(2,i) !- 0.5d0

                rr = sqrt((x-x0)**2 + (y-y0)**2)
                if(rr.lt.0.03)then
                    Us(1,i) = U00
                else
                    Us(1,i) = exp(-0.5d0*((rr-0.03)/sigma)**2) * U00
                end if
                
                Rhos(i) = RoFon * (1d0 - coef1)**(1d0/(gam-1d0))
                Us(2,i) = 0d0
                Us(3,i) = 0d0
                Ps(i) = Pfon
                Es(i) = Ps(i) / (0.4d0*Rhos(i)) + 0.5d0 * (Us(1,i)**2 + Us(2,i)**2)
                Ts(i) = TemperaturePRho(Ps(i),Rhos(i))

                Psnew(i) = Ps(i)
                Usnew(1:3,i) = Us(1:3,i)
                Esnew(i) = Es(i)
                Rhosnew(i) = Rhos(i)
                Tsnew(i) = Ts(i)
            end do
        end if

        if(InitCondInletBCon)then
            do i = 1,nCells
                r(1:3) = -meshCellCenters(1:3,i)
                call CalcNozzleBoundCond(r,U(1:3,i),T(i),P(i))

                Rho(i) = DensityPT(P(i),T(i))
                E(i) = P(i) / (0.4d0*Rho(i)) + 0.5d0 * sum(U(1:3,i)**2)

                Pnew(i) = P(i)
                Unew(1:3,i) = U(1:3,i)
                Enew(i) = E(i)
                Rhonew(i) = Rho(i)
                Tnew(i) = T(i)
            end do
            
            do i = 1,nFaces
                r(1:3) = meshFaceCenters(1:3,i)
                call CalcNozzleBoundCond(r,Us(1:3,i),Ts(i),Ps(i))

                Rhos(i) = DensityPT(Ps(i),Ts(i))
                Es(i) = Ps(i) / (0.4d0*Rhos(i)) + 0.5d0 * sum(Us(1:3,i)**2)

                Psnew(i) = Ps(i)
                Usnew(1:3,i) = Us(1:3,i)
                Esnew(i) = Es(i)
                Rhosnew(i) = Rhos(i)
                Tsnew(i) = Ts(i)
            end do
            
        end if

!	    ! conservative variables on surfaces
!	    do iCell = 1,nCells
!	        U(1,iCell) = U(1,iCell) * rho(iCell)
!	        U(2,iCell) = U(2,iCell) * rho(iCell)
!	        U(3,iCell) = U(3,iCell) * rho(iCell)
!	        E(iCell) = E(iCell) * rho(iCell)
!	    end do
    ! ************************************************************************************
	! Read Restarts
    ! ************************************************************************************
 	else ! place for restart
!        write(DirString(13:20),'(i8)') irestart
!        do i=13,20
!            if (DirString(i:i).eq.' ') DirString(i:i)='0'
!        end do
!        DirString(1:12)='Restarts/Res'

        if(autoRestart)then
            call ReadRestart()
        else
            DirString(1:3)='Res'

            if(RestartBin)then
                !call ReadRestartBin(processor(1:ie)//trim(DirString)//'/ResCFOAM.bin', &
                call ReadRestartBin(processor(1:ie)//'/ResCFOAM.bin', &
                                    ttime,nStepInitial,dt, &
                                    nCells,nFaces, &
                                    U,p,T, &
                                    Usnew,Psnew,Rhosnew, &
                                    utmp)
    !            write(*,*) '@@@ ttime,nStepInitial,dt=', ttime,nStepInitial,dt
    !            write(*,*) '@@@ nCells,nFaces=', nCells,nFaces
            else
    !            nCells = 10
    !            nFaces = 10
	            call inputScalar(processor(1:ie)//trim(DirString)//'/P',nCells,p, &
	                             ttime, nStepInitial, dt)
    !            write(*,*) ttime, nStepInitial, dt
    !            write(*,*) 'P read.', p(1),p(nCells)
	            call inputVector(processor(1:ie)//trim(DirString)//'/U',nCells,U)
    !	        write(*,*) 'U read.', U(1:3,1),U(1:3,nCells)
	            call inputVector(processor(1:ie)//trim(DirString)//'/US',nFaces,Usnew)
    !            write(*,*) 'Usnew read.'

	            call inputScalar(processor(1:ie)//trim(DirString)//'/RhoS',nFaces,Rhosnew, &
	                             ttime, nStepInitial, dt)
    !            write(*,*) 'Rhos read.'
	            call inputScalar(processor(1:ie)//trim(DirString)//'/PS',nFaces,Psnew, &
	                             ttime, nStepInitial, dt)
    !            write(*,*) 'Psnew read.'

	            call inputScalar(processor(1:ie)//trim(DirString)//'/utmp',800,utmp, &
	                             ttime, nStepInitial, dt)
    !            write(*,*) 'utmp read.'
                
	            call inputScalar(processor(1:ie)//trim(DirString)//'/T',nCells,T, &
	                             ttime, nStepInitial, dt)
    !            write(*,*) 'T read.'

    !            stop
        !              CurrentTime=ttime, nStepInitial=nStepInitial, TimeStepInitial=dt)
    	        
        !	    call inputScalar(processor(1:ie)//trim(DirString)//'/Rho',nCells,rho)
        !	    call inputScalar(processor(1:ie)//trim(DirString)//'/E',nCells,E)
        !	    call inputScalar(processor(1:ie)//trim(DirString)//'/Es',nFaces,Esnew)
        !	    call inputScalar(processor(1:ie)//trim(DirString)//'/Ts',nFaces,Tsnew)

        !	    call inputScalar(processor(1:ie)//'0/T',nCells,T)
        !	    call inputScalar(processor(1:ie)//'0/p',nCells,p)
        !	    call inputVector(processor(1:ie)//'0/U',nCells,U)
            end if
        end if

	    do iCell = 1,nCells
!	        pnew(iCell) = p(iCell)
!	        Tnew(iCell) = T(iCell)
!	        Rhonew(iCell) = Rho(iCell)
!	        Enew(iCell) = E(iCell)
!	        Unew(1,iCell) = U(1,iCell)
!	        Unew(2,iCell) = U(2,iCell)
!	        Unew(3,iCell) = U(3,iCell)
	        
!    !	p(i)=p(i)-101325.0
	        rho(iCell) = p(iCell)/(T(iCell)*cv*(Gam-1.0d0))                                     !!
!	        T(iCell) = p(iCell)/(rho(iCell)*cv*(Gam-1.0d0))                                     !!
!	        !E(iCell)= p(iCell) / ((gam-1.0d0)*rho(iCell)) + 0.5 * (U(1,iCell)**2 + U(2,iCell)**2 + U(3,iCell)**2) !/ rho(iCell)
	        !E(iCell)=p(iCell)/((gam-1.0d0)*rho(iCell))+(U(1,iCell)**2+U(2,iCell)**2+U(3,iCell)**2)/2.0d0  !!
	        E(iCell) = cv * T(iCell) +(U(1,iCell)**2+U(2,iCell)**2+U(3,iCell)**2)/2.0d0
	    end do
	    
	    do iFace=1,nFaces
	        Tsnew(iFace) = Psnew(iFace) / ((gam-1.)*cv*Rhosnew(iFace))
	        !Rhosnew(iFace) = Psnew(iFace) / ((gam-1.)*cv*Tsnew(iFace))
	        Esnew(iFace) = cv * Tsnew(iFace) + 0.5 * (Usnew(1,iFace)**2 + Usnew(2,iFace)**2 + Usnew(3,iFace)**2)
	    end do
	    
	    ! Read Pobs for Acoustic Post-Processor
	    if(AcousticsOnOff)then
   		    if (rank.eq.0) then
			    write(*,*) 'FWH: Pobs are reading ...'

                if(AcousticSurf)then
                    do iNS = 1,nNoiseSignals
                        iS = NoiseSignal(iNS)%iSurface
                        iL = NoiseSignal(iNS)%iListener
                    
                        do i = 1,300000
                            NoiseSignal(iNS)%Pobs(i) = 0.0d0
                        end do

                        nFL = NoiseSignal(iNS)%nameFileLength
                        nameFile(1:nFL-1) = NoiseSignal(iNS)%outputFile(1:nFL-1)
		                nameFile(nFL:nFL) = '1'

                        open(21,file=nameFile(1:nFL),iostat=ierr)

   		                if(ierr.ne.0) then
   		                    write(*,*) '@Warning: Acoustics Data open error!'
                            !stop 'Fatal error: Acoustics Data open error!'
                        else
                            do i=1,300000
			                    read(21,*,iostat=ierr) auxReadReal,NoiseSignal(iNS)%Pobs(i)
		                        if(ierr.lt.0)then ! EOF
			                        write(*,*) 'FWH: Pobs was read ', i, '!'
			                        exit
                                end if
                                if(ierr.gt.0)then ! error during the read execution
                                    stop '@@Fatal error: Acoustics Data is incorrect!'            
                                end if
		                    end do
		                end if
		                
		                close(21)
		            end do
                else
   		            open(20,file='Pobs.csv',iostat=ierr)
   		            if(ierr.ne.0) then
                        stop 'Fatal error: Pobs.csv open error!'
                    end if
                    do i = 1,300000
                        Pobs(i) = 0.0d0
                    end do

			        do i=1,300000
				        read(20,*,iostat=ierr) auxReadReal,Pobs(i)
		                if(ierr.lt.0)then ! EOF
			                write(*,*) '@Warning: FWH: Pobs was not read full!'
			                exit
                        end if
                        if(ierr.gt.0)then ! error during the read execution
                            stop '@@Fatal error: Pobs.csv is incorrect!'            
                        end if
			        end do
			        close(20)
			    end if
			    write(*,*) 'FWH: Pobs was read successfully'
		    end if
        end if

!        do iTrans = 1,nTransducers
!            Transducer(iTrans)%TimeStart = max(Transducer(iTrans)%TimeStart,ttime)
!        end do

	end if
    ! ************************************************************************************

    dt = TimeStepMin ! min(TimeStepMax,max(dt,TimeStepMin))
	dt2=dt/2.0d0

	if(AsynchronousTime)then
	    allocate(SubStepCellVars(5,nCells))
	    allocate(CharsCellPresent(nCells))
	    
!	    allocate(CharsFacePresent(2,nFaces))

        allocate(CharsOld(5,2,nFaces),CharsNew(5,2,nFaces))
        allocate(fluxFaces(5,2,nFaces))
        allocate(fluxFacesReserve(5,2,nFaces))
        allocate(fluxesFacePresent(2,nFaces))
        allocate(FaceVarsNew(5,2,nFaces), FaceVarsOld(5,2,nFaces))
        
        allocate(fluxCorrection(5,2,nFaces))

        allocate(fluxesCellPresent(nCells))
        allocate(sumCellFluxes(5,nCells))

        if(CalcViscosityOn)then
            allocate(CurrentVelocities(3,nCells))
        end if
        allocate(ViscousFluxes(3,2,nFAces))

        allocate(CellGlobalStepDone(nCells))

        allocate(TimeCell(nCells),dtCell(nCells),dt2Cell(nCells),NextTimeCell(nCells))
        allocate(nNextTimeCell(nCells),nNextTimeCell1(nCells))
        allocate(dnTimeCell(nCells))
!        allocate(CurrentTimeCell(nCells))
        allocate(dtFace(nFaces))
!        allocate(dtFaceCorr(2,nFaces))

        allocate(Acoef(2,nFaces))
        allocate(AcoefCell(nCells))

        allocate(CellViscousStencil(nCells))
        allocate(FaceViscousStencil(nFaces))
        allocate(PointViscousStencil(nPoints))

!        CharsCellPresent = .true.
        fluxesFacePresent = .false.
        
!        CharsFacePresent = .false.

        CellGlobalStepDone = .false.

        SubStepCellVars = 0d0
        fluxFacesReserve = 0d0
        fluxFaces = 0d0
        CharsOld = 0d0
        CharsNew = 0d0
        sumCellFluxes = 0d0
        fluxCorrection = 0d0
        ViscousFluxes = 0d0

        TimeLocal = 0d0 ! ttime
        TimeCell(1:nCells) = 0d0 ! ttime
        
        nTimeStepDiff = int(TimeStepDiff/TimeStepMin)
!        CurrentTimeCell(1:nCells) = ttime

        dtCell(1:nCells) = dt ! TimeStepInitial
        dt2Cell(1:nCells) = 0.5d0 * dtCell(1:nCells)
        
!        NextTimeCell(1:nCells) = ttime + dtCell(1:nCells)
!        nNextTimeCell(1:nCells) = int((ttime + dtCell(1:nCells))/TimeStepMin)
!        NextTimeLocal = minVal(NextTimeCell(1:nCells))
!        nNextTimeLocal = minVal(nNextTimeCell(1:nCells))

!! Calculation InitCond through Calc3D
        if(InitCondCalc3D) then
            do iCell = 1,nCells
                x = meshCellCenters(1,iCell)
                c = 340d0
                TimeX = x / c
                r = 0d0
                call CalcQ(r,TimeX,Q,dQdt,p0,U0)

                r(1:3) = -meshCellCenters(1:3,iCell)
                rMod = sqrt(sum(r(1:3)*r(1:3)))
                c = sqrt(1.4*PFon/RoFon)
                TimeX = - rMod / c
                call CalcQ3D(r,0d0,Q,dQdt,p0,U0)

                Pnew(iCell) = p0
                Unew(1:3,iCell) = U0(1:3)
                Rhonew(iCell) = RoFon + (P0-PFon)/c**2
        	    Enew(iCell) = P0 / (0.4 * (RoFon + (P0-PFon)/c**2)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
    !            Pnew(iCell) = 0.4d0 * (Enew(iCell) - 0.5d0 * sum(U0(1:3)*U0(1:3))) * Rhonew(iCell)! * 0.4d0

                p(iCell) = Pnew(iCell)
                U(1:3,iCell) = Unew(1:3,iCell)
                E(iCell) = Enew(iCell)
                Rho(iCell) = Rhonew(iCell)
                
    !            if(r(1).gt.0d0.and.r(2).gt.0d0.and.r(3).lt.0d0) then
    !                write(771,'(8f20.10)') r(1:3), sqrt(sum(r(1:3)*r(1:3))), p0, U0(1:3)
    !            end if
                
            end do

            do iFace = 1,nFaces
                x = -meshFaceCenters(1,iFace)
                c = 340d0
                TimeX = x / c
                call CalcQ(r,TimeX,Q,dQdt,p0,U0)

                r(1:3) = meshFaceCenters(1:3,iFace)
                rMod = sqrt(sum(r(1:3)*r(1:3)))
                c = sqrt(1.4*PFon/RoFon)
                TimeX = - rMod / c
                call CalcQ3D(r,0d0,Q,dQdt,p0,U0)

                Psnew(iFace) = p0
                Usnew(1:3,iFace) = U0(1:3)
        	    Esnew(iFace) = P0 / (0.4 * (RoFon + (P0-PFon)/c**2)) + 0.5d0 * sum(U0(1:3)*U0(1:3))
                Rhosnew(iFace) = RoFon + (P0-PFon)/c**2
    !            Psnew(iFace) = (Esnew(iFace) - 0.5d0 * sum(U0(1:3)*U0(1:3))) * Rhosnew(iFace) * 0.4d0

                ps(iFace) = Psnew(iFace)
                Us(1:3,iFace) = Usnew(1:3,iFace)
                Es(iFace) = Esnew(iFace)
                Rhos(iFace) = Rhosnew(iFace)

    !            if(r(1).gt.0d0.and.r(2).gt.0d0.and.r(3).lt.0d0) then
    !                write(771,'(8f20.10)') r(1:3), sqrt(sum(r(1:3)*r(1:3))), p0, U0(1:3)
    !            end if

    !            if(abs(r(2)).lt.1d1.and.abs(r(3)).lt.1d1) then
    !                write(771,'(8f20.10)') r(1:3), sqrt(sum(r(1:3)*r(1:3))), p0, U0(1:3)
    !            end if
            end do
        end if
!-------------------------------------------------------------------!

!!! Calculation InitCond through BlowoffInitialCondition(r,p,U,Rho)
!        do iCell = 1,nCells
!
!            r(1:3) = -meshCellCenters(1:3,iCell)
!            call BlowoffInitialCondition(r,Pnew(iCell),Unew(1:3,iCell),Rhonew(iCell))
!            
!        	Enew(iCell) = EnergyPRhoU(Pnew(iCell),Rhonew(iCell),Unew(1:3,iCell))
!
!            p(iCell) = Pnew(iCell)
!            U(1:3,iCell) = Unew(1:3,iCell)
!            E(iCell) = Enew(iCell)
!            Rho(iCell) = Rhonew(iCell)
!        end do
!
!        do iFace = 1,nFaces
!            r(1:3) = meshFaceCenters(1:3,iFace)
!            call BlowoffInitialCondition(r,Psnew(iCell),Usnew(1:3,iCell),Rhosnew(iCell))
!            
!        	Esnew(iCell) = EnergyPRhoU(Pnew(iCell),Rhonew(iCell),Unew(1:3,iCell))
!
!            ps(iFace) = Psnew(iFace)
!            Us(1:3,iFace) = Usnew(1:3,iFace)
!            Es(iFace) = Esnew(iFace)
!            Rhos(iFace) = Rhosnew(iFace)
!        end do
!!-------------------------------------------------------------------!
        do iCell = 1,nCells
            SubStepCellVars(1,iCell) = rho(iCell)
            SubStepCellVars(2:4,iCell) = U(1:3,iCell)
            SubStepCellVars(5,iCell) = E(iCell)
        end do

        do iFace = 1,nFaces
            do i = 1,2
                FaceVarsOld(1,i,iFace) = Rhosnew(iFace)
                do j = 1,3
                    FaceVarsOld(j+1,i,iFace) = Usnew(j,iFace)
!                    FaceVarsOld(2:4,2,iFace) = Usnew(1:3,iFace)
                end do
                FaceVarsOld(5,i,iFace) = Psnew(iFace) ! Esnew(iFace)
            end do
        end do

        do iFace = 1,nFaces
            do i = 1,2
                do j = 1,5
                    FaceVarsNew(j,i,iFace) = FaceVarsOld(j,i,iFace)
!        FaceVarsNew(1:5,1:2,1:nFaces) = FaceVarsOld(1:5,1:2,1:nFaces)
                end do
            end do
        end do

        do iCell = 1,nCells
            RhoNew(iCell) = Rho(iCell)
            do i = 1,3
                Unew(i,iCell) = U(i,iCell)
            end do
            Enew(iCell) = E(iCell)
            Pnew(iCell) = P(iCell)
        end do

	    if(CalcViscosityOn)then
	        do iCell = 1,nCells
	            do i =1,3
	                CurrentVelocities(i,iCell) = Unew(i,iCell)
	            end do
	        end do
	    end if

        FaceViscousStencil = .true.
        PointViscousStencil = .true.
        nNextTimeCell1 = .false.

        call CalcCharsInit()
        call CalcFluxesInit()
        call CalcSumFluxes()
        
!        fluxesFacePresent = .true.
!        CharsCellPresent(iCell) = .false.
!        call CalcSubStepVars()



!	    call OutputTransducerAsynData()

	end if


!***************************************************************************
    ! end set initial conditions
!***************************************************************************

    if(rank.eq.0)then
        write(*,*) 'Initialization has been done successfully!'
    end if

end subroutine SetInitConditions

subroutine CalcCharsInit()
use Data_Mesh, only: nFaces, meshFaceOwner,meshFaceNeighbour, meshSf
use DataFaces, only: CharsNew,CharsOld
use DataCells, only: SubStepCellVars
use DataPerfectGas, only: PressureAllVars, SoundVelocityPRho, ACoefPRho
implicit none

    integer :: iFace, iCell
    real(8) :: EigenVal3,Sound, Pres,Rho, Acoef
    real(8) :: a(3,3), n(3),normal(3)
    
    do iFace = 1,nFaces
        ! Definition of Characteristics, Eigen Values and Sound Velocities in
        n(1:3) = meshSf(1:3,iFace)
        normal(1:3) = n(1:3)/dsqrt(sum(n(1:3)*n(1:3)))

        ! transformation {x,y,z} -> {normal,tan1,tan2}
        call coor(normal(1),normal(2),normal(3),a)

        iCell = meshFaceOwner(iFace)
        Pres = PressureAllVars(SubStepCellVars(1:5,iCell))
        Rho = SubStepCellVars(1,iCell)
        Acoef = ACoefPRho(Pres,Rho) ! GammY * Pres * Rho**(-Gamma)
        call FaceVars2Chars(iFace,iFace,a,Acoef, EigenVal3,Sound,1, CharsOld(1:5,1,iFace))
        
        iCell = meshFaceOwner(iFace)
        Pres = PressureAllVars(SubStepCellVars(1:5,iCell))
        Rho = SubStepCellVars(1,iCell)
        Acoef = ACoefPRho(Pres,Rho) ! GammY * Pres * Rho**(-Gamma)
        call FaceVars2Chars(iFace,iFace,a,Acoef, EigenVal3,Sound,1, CharsOld(1:5,2,iFace))
    end do

end subroutine CalcCharsInit

subroutine CalcFluxesInit()
use Data_Cfoam, only: Usnew, Esnew, Psnew, Rhosnew
use Data_Mesh, only: nFaces, meshSf
use DataFaces,only: fluxesFacePresent, fluxFaces
implicit none

    integer :: iFace
    real(8) :: vtmp
    real(8) :: Fluxes(5)

    do iFace = 1,nFaces
        if(fluxesFacePresent(1,iFace).and.fluxesFacePresent(2,iFace)) cycle
        
        vtmp = sum(meshSf(1:3,iFace)*Usnew(1:3,iFace))
        Fluxes(1) = vtmp * Rhosnew(iFace)
        Fluxes(2:4) = vtmp*Rhosnew(iFace)*Usnew(1:3,iFace) + (Psnew(iFace)-1d5)*meshSf(1:3,iFace)
!        Fluxes(2:4) = vtmp*Rhosnew(iFace)*Usnew(1:3,iFace) + Psnew(iFace)*meshSf(1:3,iFace)
        Fluxes(5) = vtmp * (Rhosnew(iFace)*Esnew(iFace) + Psnew(iFace))

        if(fluxesFacePresent(1,iFace).eqv..false.)then
            fluxFaces(1:5,1,iFace) = Fluxes(1:5)
            fluxesFacePresent(1,iFace) = .true.
        end if
        if(fluxesFacePresent(2,iFace).eqv..false.)then
            fluxFaces(1:5,2,iFace) = -Fluxes(1:5)
            fluxesFacePresent(2,iFace) = .true.
        end if
    end do

end subroutine CalcFluxesInit

subroutine CalcSumFluxes()
use Data_Mesh, only: nCells, meshCells, meshFaceOwner
use DataFaces, only: fluxesFacePresent, fluxFaces, ViscousFluxes
use DataCells, only: sumCellFluxes, fluxesCellPresent
implicit none

    integer :: iCell, i, numFace, iFace
    logical :: CellFluxesPresent
    real(8) :: Fluxes(5)

    do iCell = 1,nCells
        fluxesCellPresent(iCell) = .true.
        Fluxes(1:5) = 0d0
        do numFace = 1,6
            iFace = meshCells(numFace,iCell)
            if(meshFaceOwner(iFace).eq.iCell)then
                i = 1
            else
                i = 2
            end if
            if(fluxesFacePresent(i,iFace))then
                Fluxes(1:5) = Fluxes(1:5) + fluxFaces(1:5,i,iFace)
!                Fluxes(2:4) = Fluxes(2:4) + ViscousFluxes(1:3,i,iFace)
            else
                fluxesCellPresent(iCell) = .false.
                exit
            end if
        end do
        if(fluxesCellPresent(iCell))then
            sumCellFluxes(1:5,iCell) = Fluxes(1:5)
        end if
    end do

end subroutine CalcSumFluxes

!subroutine CalcSubStepVars()
!use Data_Mesh, only: nCells, meshV
!use DataCells, only: sumCellFluxes, fluxesCellPresent, SubStepCellVars, CharsCellPresent
!use DataTime, only: dt2Cell
!implicit none
!
!    integer :: iCell
!
!    do iCell = 1,nCells
!!        if(fluxesCellPresent(iCell).eq..false.) cycle
!
!        SubStepCellVars(2:5,iCell) = SubStepCellVars(1,iCell) * SubStepCellVars(2:5,iCell) - &
!                                     sumCellFluxes(2:5,iCell) * dt2Cell(iCell) / meshV(iCell)
!        SubStepCellVars(1,iCell) = SubStepCellVars(1,iCell) - &
!                                     sumCellFluxes(1,iCell) * dt2Cell(iCell) / meshV(iCell)
!        SubStepCellVars(2:5,iCell) = SubStepCellVars(2:5,iCell) / SubStepCellVars(1,iCell)
!    end do
!
!end subroutine CalcSubStepVars
